## The basic files and libraries needed for most presentations
# creates the libraries and common-functions sections
read_chunk("../common/utility_functions.R")

require(ggplot2) #for plots
require(lattice) # nicer scatter plots
require(plyr) # for processing data.frames
require(grid) # contains the arrow function
require(biOps) # for basic image processing
require(doMC) # for parallel code
require(png) # for reading png images
require(gridExtra)
require(reshape2) # for the melt function
## To install EBImage
# source("http://bioconductor.org/biocLite.R")
# biocLite("EBImage")
require(EBImage) # for more image processing
used.libraries<-c("ggplot2","lattice","plyr","reshape2","grid","gridExtra","biOps","png","EBImage")

# start parallel environment
registerDoMC()
# functions for converting images back and forth
im.to.df<-function(in.img) {
  out.im<-expand.grid(x=1:nrow(in.img),y=1:ncol(in.img))
  out.im$val<-as.vector(in.img)
  out.im
}
df.to.im<-function(in.df,val.col="val",inv=F) {
  in.vals<-in.df[[val.col]]
  if(class(in.vals[1])=="logical") in.vals<-as.integer(in.vals*255)
  if(inv) in.vals<-255-in.vals
  out.mat<-matrix(in.vals,nrow=length(unique(in.df$x)),byrow=F)
  attr(out.mat,"type")<-"grey"
  out.mat
}
ddply.cutcols<-function(...,cols=1) {
  # run standard ddply command 
  cur.table<-ddply(...)
  cutlabel.fixer<-function(oVal) {
    sapply(oVal,function(x) {
      cnv<-as.character(x)
      mean(as.numeric(strsplit(substr(cnv,2,nchar(cnv)-1),",")[[1]]))
    })
  }
  cutname.fixer<-function(c.str) {
    s.str<-strsplit(c.str,"(",fixed=T)[[1]]
    t.str<-strsplit(paste(s.str[c(2:length(s.str))],collapse="("),",")[[1]]
    paste(t.str[c(1:length(t.str)-1)],collapse=",")
  }
  for(i in c(1:cols)) {
    cur.table[,i]<-cutlabel.fixer(cur.table[,i])
    names(cur.table)[i]<-cutname.fixer(names(cur.table)[i])
  }
  cur.table
}

show.pngs.as.grid<-function(file.list,title.fun,zoom=1) {
  preparePng<-function(x) rasterGrob(readPNG(x,native=T,info=T),width=unit(zoom,"npc"),interp=F)
  labelPng<-function(x,title="junk") (qplot(1:300, 1:300, geom="blank",xlab=NULL,ylab=NULL,asp=1)+
                                        annotation_custom(preparePng(x))+
                                        labs(title=title)+theme_bw(24)+
                                        theme(axis.text.x = element_blank(),
                                              axis.text.y = element_blank()))
  imgList<-llply(file.list,function(x) labelPng(x,title.fun(x)) )
  do.call(grid.arrange,imgList)
}
## Standard image processing tools which I use for visualizing the examples in the script
commean.fun<-function(in.df) {
  ddply(in.df,.(val), function(c.cell) {
    weight.sum<-sum(c.cell$weight)
    data.frame(xv=mean(c.cell$x),
               yv=mean(c.cell$y),
               xm=with(c.cell,sum(x*weight)/weight.sum),
               ym=with(c.cell,sum(y*weight)/weight.sum)
    )
  })
}

colMeans.df<-function(x,...) as.data.frame(t(colMeans(x,...)))

pca.fun<-function(in.df) {
  ddply(in.df,.(val), function(c.cell) {
    c.cell.cov<-cov(c.cell[,c("x","y")])
    c.cell.eigen<-eigen(c.cell.cov)
    
    c.cell.mean<-colMeans.df(c.cell[,c("x","y")])
    out.df<-cbind(c.cell.mean,
                  data.frame(vx=c.cell.eigen$vectors[1,],
                             vy=c.cell.eigen$vectors[2,],
                             vw=sqrt(c.cell.eigen$values),
                             th.off=atan2(c.cell.eigen$vectors[2,],c.cell.eigen$vectors[1,]))
    )
  })
}
vec.to.ellipse<-function(pca.df) {
  ddply(pca.df,.(val),function(cur.pca) {
    # assume there are two vectors now
    create.ellipse.points(x.off=cur.pca[1,"x"],y.off=cur.pca[1,"y"],
                          b=sqrt(5)*cur.pca[1,"vw"],a=sqrt(5)*cur.pca[2,"vw"],
                          th.off=pi/2-atan2(cur.pca[1,"vy"],cur.pca[1,"vx"]),
                          x.cent=cur.pca[1,"x"],y.cent=cur.pca[1,"y"])
  })
}

# test function for ellipse generation
# ggplot(ldply(seq(-pi,pi,length.out=100),function(th) create.ellipse.points(a=1,b=2,th.off=th,th.val=th)),aes(x=x,y=y))+geom_path()+facet_wrap(~th.val)+coord_equal()
create.ellipse.points<-function(x.off=0,y.off=0,a=1,b=NULL,th.off=0,th.max=2*pi,pts=36,...) {
  if (is.null(b)) b<-a
  th<-seq(0,th.max,length.out=pts)
  data.frame(x=a*cos(th.off)*cos(th)+b*sin(th.off)*sin(th)+x.off,
             y=-1*a*sin(th.off)*cos(th)+b*cos(th.off)*sin(th)+y.off,
             id=as.factor(paste(x.off,y.off,a,b,th.off,pts,sep=":")),...)
}
deform.ellipse.draw<-function(c.box) {
  create.ellipse.points(x.off=c.box$x[1],
                        y.off=c.box$y[1],
                        a=c.box$a[1],
                        b=c.box$b[1],
                        th.off=c.box$th[1],
                        col=c.box$col[1])                    
}
bbox.fun<-function(in.df) {
  ddply(in.df,.(val), function(c.cell) {
    c.cell.mean<-colMeans.df(c.cell[,c("x","y")])
    xmn<-emin(c.cell$x)
    xmx<-emax(c.cell$x)
    ymn<-emin(c.cell$y)
    ymx<-emax(c.cell$y)
    out.df<-cbind(c.cell.mean,
                  data.frame(xi=c(xmn,xmn,xmx,xmx,xmn),
                             yi=c(ymn,ymx,ymx,ymn,ymn),
                             xw=xmx-xmn,
                             yw=ymx-ymn
                  ))
  })
}

# since the edge of the pixel is 0.5 away from the middle of the pixel
emin<-function(...) min(...)-0.5
emax<-function(...) max(...)+0.5
extents.fun<-function(in.df) {
  ddply(in.df,.(val), function(c.cell) {
    c.cell.mean<-colMeans.df(c.cell[,c("x","y")])
    out.df<-cbind(c.cell.mean,data.frame(xmin=c(c.cell.mean$x,emin(c.cell$x)),
                                         xmax=c(c.cell.mean$x,emax(c.cell$x)),
                                         ymin=c(emin(c.cell$y),c.cell.mean$y),
                                         ymax=c(emax(c.cell$y),c.cell.mean$y)))
  })
}

common.image.path<-"../common"
qbi.file<-function(file.name) file.path(common.image.path,"figures",file.name)
qbi.data<-function(file.name) file.path(common.image.path,"data",file.name)

th_fillmap.fn<-function(max.val) scale_fill_gradientn(colours=rainbow(10),limits=c(0,max.val))

grad.beam.profile<-matrix(rep(0.9+0.00*c(1:nrow(gray.img)),ncol(gray.img)),dim(gray.img))
sample.img<-1.0*exp(-material.img*0.1)
ccd.img<-sample.img * grad.beam.profile
model.image<-subset(
  rbind(
  cbind(im.to.df(grad.beam.profile),type="Beam"),
  cbind(im.to.df(sample.img),type="Sample"),
  cbind(im.to.df(ccd.img),type="Detector")
  ),
  x%%3==0 & y%%3==0 # downsample by 9
  )
ggplot(model.image,aes(x=y,y=x,fill=val))+
  geom_tile()+
  guides(fill=F)+
  facet_grid(type~.)+
  coord_equal()+
  theme_bw(20)

grad.beam.profile<-matrix(rep(0.9+0.001*c(1:nrow(gray.img)),ncol(gray.img)),dim(gray.img))
sample.img<-1.0*exp(-material.img*0.1)
ccd.img<-sample.img * grad.beam.profile
model.image<-subset(
  rbind(
  cbind(im.to.df(grad.beam.profile),type="Beam"),
  cbind(im.to.df(sample.img),type="Sample"),
  cbind(im.to.df(ccd.img),type="Detector")
  ),
  x%%3==0 & y%%3==0 # downsample by 9
  )
ggplot(model.image,aes(x=y,y=x,fill=val))+
  geom_tile()+
  facet_grid(type~.)+
  coord_equal()+
  labs(fill="")+
  scale_fill_gradientn(colours=rainbow(6))+
  theme_bw(20)

beam.profile<-mutate(
  expand.grid(x=1:nrow(material.img),y=1:ncol(material.img)),
  val=exp(-((x-mean(x))^2+(y-mean(y))^2)/350^2) # just a simple gaussian profile
  )
beam.profile<-df.to.im(beam.profile)
sample.img<-1.0*exp(-material.img*0.1)
ccd.img<-sample.img * beam.profile
model.image<-subset(
  rbind(
  cbind(im.to.df(beam.profile),type="Beam"),
  cbind(im.to.df(sample.img),type="Sample"),
  cbind(im.to.df(ccd.img),type="Detector")
  ),
  x%%3==0 & y%%3==0 # downsample by 9
  )
ggplot(model.image,aes(x=y,y=x,fill=val))+
  geom_tile()+
  facet_grid(type~.)+
  coord_equal()+
  scale_fill_gradientn(colours=rainbow(6))+
  labs(fill="")+
  theme_bw(20)

nx<-4
ny<-4
grad.im<-expand.grid(x=c(-nx:nx)/nx*2*pi,y=c(-ny:ny)/ny*2*pi)
phase.im<-runif(nrow(grad.im))
grad.im<-cbind(grad.im,
               col=1*(phase.im>0.66)+
                 2*(phase.im<0.33)+
                 0.4*runif(nrow(grad.im),min=-1,max=1))
bn.wid<-diff(range(grad.im$col))/10

bl.fun<-function(x) 1*exp(-x*0.5)

ggplot(grad.im,aes(x=col,y=bl.fun(col)))+
  geom_point(aes(color=cut_interval(col,3)))+
  geom_segment(data=data.frame(
    x=c(0.5,1.5),y=c(0.25),
    xend=c(0.5,1.5),yend=c(1.25)
    ),aes(x=x,y=y,xend=xend,yend=yend),fill="red",size=3,alpha=0.5)+
  labs(x=expression(alpha(x,y)),y=expression(paste("I"[sample],"(x,y)")),title="Probability Density Function",fill="Groups")+theme_bw(20)

ggplot(grad.im,aes(x=col))+
  geom_density(aes(fill=cut_interval(col,3)))+
  geom_segment(data=data.frame(
    x=c(0.5,1.5),y=c(0),
    xend=c(0.5,1.5),yend=c(6)
    ),aes(x=x,y=y,xend=xend,yend=yend),fill="red",size=3,alpha=0.5)+
  labs(x=expression(alpha(x,y)),y="Frequency",title="Probability Density Function",fill="Groups")+theme_bw(20)

nx<-5
ny<-5
cross.im<-expand.grid(x=c(-nx:nx)/nx*2*pi,y=c(-ny:ny)/ny*2*pi)
cross.im<-cbind(cross.im,
               col=1.5*with(cross.im,abs(cos(x*y))/(abs(x*y)+(3*pi/nx)))+
                 0.5*runif(nrow(cross.im)))
bn.wid<-diff(range(cross.im$col))/10
ggplot(cross.im,aes(x=x,y=y,fill=col))+
  geom_tile()+
  scale_fill_gradient(low="black",high="white")+
  labs(fill="Intensity")+
  theme_bw(20)

ggplot(cross.im,aes(x=col))+geom_histogram(binwidth=bn.wid)+
  labs(x="Intensity",title="Probability Density Function")+
  theme_bw(20)

thresh.val<-0.75
cross.im$val<-(cross.im$col>=thresh.val)
ggplot(cross.im,aes(x=x,y=y))+
  geom_tile(aes(fill=col))+
  geom_tile(data=subset(cross.im,val),fill="red",color="black",alpha=0.3)+
  geom_tile(data=subset(cross.im,!val),fill="blue",color="black",alpha=0.3)+
  scale_fill_gradient(low="black",high="white")+
  labs(fill="Intensity")+
  theme_bw(20)

ggplot(cbind(cross.im,in.thresh=cross.im$col>=thresh.val),aes(x=col))+
  geom_histogram(binwidth=bn.wid,aes(fill=in.thresh))+
  labs(x="Intensity",color="In Threshold")+scale_fill_manual(values=c("blue","red"))+
  theme_bw(20)

thresh.vals<-c(2:10)/10
grad.im.th<-ldply(thresh.vals,function(thresh.val) 
  cbind(cross.im,thresh=thresh.val,in.thresh=(cross.im$col>=thresh.val)))
ggplot(grad.im.th,aes(x=col))+
  geom_histogram(binwidth=bn.wid,aes(fill=in.thresh))+
  labs(x="Intensity",fill="Above Threshold")+facet_wrap(~thresh)+
  theme_bw(15)

ggplot(grad.im.th,aes(x=x,y=y))+
  geom_tile(aes(fill=in.thresh),color="black",alpha=0.75)+
  labs(fill="Above Threshold")+facet_wrap(~thresh)+
  theme_bw(20)

cellImage<-im.to.df(jpeg::readJPEG("ext-figures/Cell_Colony.jpg"))
ggplot(cellImage,aes(x=x,y=y,fill=val))+
  geom_tile()+
  labs(fill="Intensity")+coord_equal()+
  theme_bw(20)

th.vals<-seq(0.4,0.85,length.out=3)
thlabel<-function(x,...) switch(x,"Too low","Good","Too high")
im.vals<-ldply(1:length(th.vals),function(th.val)
  cbind(cellImage,
        in.thresh=ifelse(cellImage$val

nx<-7
ny<-7
n.pi<-4
grad.im<-expand.grid(x=c(-nx:nx)/nx*n.pi*pi,y=c(-ny:ny)/ny*n.pi*pi)

grad.im<-cbind(grad.im,
               col=1.5*with(grad.im,abs(cos(x*y))/(abs(x*y)+(3*pi/nx)))+
                 0.5*runif(nrow(grad.im)),
               x.vec=with(grad.im,y),
               y.vec=with(grad.im,x))
# normalize vector
grad.im[,c("x.vec","y.vec")]<-with(grad.im,cbind(x.vec/(sqrt(x.vec^2+y.vec^2)),
                                         y.vec/(sqrt(x.vec^2+y.vec^2))))
bn.wid<-c(diff(range(grad.im$x.vec))/10,diff(range(grad.im$y.vec))/10)
ggplot(grad.im,aes(x=x,y=y,fill=col))+
  geom_segment(aes(xend=x+x.vec,yend=y+y.vec),arrow=arrow(length = unit(0.15,"cm")),size=2)+
  scale_fill_gradient(low="black",high="white")+
  guides(fill=F)+labs(title="Orientation Map")+
  theme_bw(20)

ggplot(grad.im,aes(x=x.vec,y=y.vec))+
  stat_bin2d(binwidth=c(0.25,.25),drop=F)+
  labs(x="Pfx",y="Pfy",fill="Frequency",title="Orientation Histogram")+
  xlim(-1,1)+ylim(-1,1)+
  theme_bw(20)

g.with.thresh<-cbind(grad.im,in.thresh=with(grad.im,x.vec>0.25 & y.vec>0.25))
ggplot(g.with.thresh,
       aes(x=x,y=y,fill=in.thresh,color=in.thresh))+
  geom_tile(alpha=0.5,aes(fill=in.thresh))+
  geom_segment(aes(xend=x+x.vec,yend=y+y.vec),arrow=arrow(length = unit(0.2,"cm")),size=3)+
  labs(color="In Threshold")+guides(fill=FALSE)+
  theme_bw(20)

bn.wid<-c(0.25,.25)
keep.bins<-expand.grid(x.vec=seq(-1,1,bn.wid[1]/10),y.vec=seq(-1,1,bn.wid[2]/10))
keep.bins<-cbind(keep.bins,in.thresh=with(keep.bins,x.vec>0.25 & y.vec>0.25))
ggplot(g.with.thresh,aes(x=x.vec,y=y.vec))+
  stat_bin2d(binwidth=bn.wid,drop=F)+
  geom_tile(data=subset(g.with.thresh,in.thresh),fill="red",alpha=0.4)+
  labs(x="Pfx",y="Pfy",fill="Threshold")+xlim(-1,1)+ylim(-1,1)+
  theme_bw(20)

i.vec<-c(1,0)
j.vec<-c(0,1)
g.cmp.thresh<-cbind(grad.im,in.thresh=with(grad.im,
                                            abs(x.vec*i.vec[1]+y.vec*i.vec[2])==1
                                            ))
ggplot(g.cmp.thresh,
       aes(x=x,y=y,fill=in.thresh,color=in.thresh))+
  geom_tile(alpha=0.5,aes(fill=in.thresh))+
  geom_segment(aes(xend=x+x.vec,yend=y+y.vec),arrow=arrow(length = unit(0.2,"cm")),size=3)+
  labs(color="In Threshold")+guides(fill=FALSE)+
  theme_bw(20)

i.vec<-c(1,0)
j.vec<-c(0,1)
ang.accept<-function(c.ang) g.cmp.thresh<-cbind(grad.im,
                                                ang.val=c.ang,
                                                in.thresh=with(grad.im,
                                            acos(x.vec*i.vec[1]+y.vec*i.vec[2])<=c.ang/180*pi
                                            ))
ang.vals<-seq(5,180,length.out=9)
ggplot(ldply(ang.vals,ang.accept),
       aes(x=x,y=y,fill=in.thresh,color=in.thresh))+
  geom_tile(alpha=0.5,aes(fill=in.thresh))+
  geom_segment(aes(xend=x+x.vec,yend=y+y.vec),arrow=arrow(length = unit(0.2,"cm")),size=3)+
  facet_wrap(~ang.val)+
  labs(color="In Threshold")+guides(fill=FALSE)

  theme_bw(20)

nx<-4
ny<-4
n.pi<-4
grad.im<-expand.grid(x=c(-nx:nx)/nx*n.pi*pi,
                     y=c(-ny:ny)/ny*n.pi*pi)

grad.im<-cbind(grad.im,
               col=1.5*with(grad.im,abs(cos(x*y))/(abs(x*y)+(3*pi/nx)))+
                 0.5*runif(nrow(grad.im)),
              a=with(grad.im,sqrt(2/(abs(x)+0.5))),
               b=with(grad.im,0.5*sqrt(abs(x)+1)),
              th=0.5*runif(nrow(grad.im)),
              aiso=1,count=1)

create.ellipse.points<-function(x.off=0,y.off=0,a=1,b=NULL,th.off=0,th.max=2*pi,pts=36,...) {
  if (is.null(b)) b<-a
  th<-seq(0,th.max,length.out=pts)
  data.frame(x=a*cos(th.off)*cos(th)+b*sin(th.off)*sin(th)+x.off,
             y=-1*a*sin(th.off)*cos(th)+b*cos(th.off)*sin(th)+y.off,
             id=as.factor(paste(x.off,y.off,a,b,th.off,pts,sep=":")),...)
}
deform.ellipse.draw<-function(c.box) {
  create.ellipse.points(x.off=c.box$x[1],
                        y.off=c.box$y[1],
                        a=c.box$a[1],
                        b=c.box$b[1],
                        th.off=c.box$th[1],
                        col=c.box$col[1])                    
}

# normalize vector
tens.im<-ddply(grad.im,.(x,y),deform.ellipse.draw)

ggplot(tens.im,aes(x=x,y=y,group=as.factor(id),fill=col))+
  geom_polygon(color="black")+coord_fixed(ratio=1)+scale_fill_gradient(low="black",high="white")+guides(fill=F)+
  theme_bw(20)

ggplot(grad.im)+
  geom_histogram(aes(x=a,fill="Width"),alpha=0.7)+
  geom_histogram(aes(x=b,fill="Height"),alpha=0.7)+
  geom_histogram(aes(x=th,fill="Orientation"),alpha=0.7)+
  geom_histogram(aes(x=col,fill="Color"),alpha=0.7)+
  guides(color=F)+labs(fill="Variable")+
  theme_bw(15)

plot(grad.im[,c("a","b","th","col")])

ggplot(shale.vals,aes(x=x,y=y,fill=th.label,alpha=1-val))+coord_fixed(ratio=1)+
  geom_raster()+guides(fill=F,alpha=F)+theme_bw(20)+facet_wrap(~th.label)

make.circle<-function(nx,r,nscale=1,
                      circ.fun=function(x,y,r) ((x^2+y^2)

ggplot(sph.list,aes(x=val,color=as.factor(r)))+
  geom_density(trim=T,adjust=1,kernel="biweight")+#facet_wrap(~r)+#scale_y_log10()+
  theme_bw(20)+labs(x="Intensity",y="Frequency",color="Radius")

nx<-8
ny<-8
cross.im<-expand.grid(x=c(-nx:nx)/nx*2*pi,y=c(-ny:ny)/ny*2*pi)
cross.im<-cbind(cross.im,
               col=with(cross.im,1.5*((abs(x)<2) | (abs(y)<2))+
                 0.5*runif(nrow(cross.im))))
thresh.val<-0.75
cross.im$val<-(cross.im$col>=thresh.val)
ggplot(cross.im,aes(x=x,y=y))+
  geom_tile(aes(fill=col))+
  geom_tile(data=subset(cross.im,val),fill="red",color="black",alpha=0.3)+
  geom_tile(data=subset(cross.im,!val),fill="blue",color="black",alpha=0.3)+
  scale_fill_gradient(low="black",high="white")+
  labs(fill="Intensity")+
  theme_bw(20)

morph.2d<-expand.grid(x=c(-1,0,1),y=c(-1,0,1))
morph.2d$r<-with(morph.2d,sqrt(x^2+y^2))

ggplot(morph.2d,aes(x=x,y=y))+
  geom_tile(color="black",
            aes(fill=ifelse(r==0,"Pixel of Interest",
                                   ifelse(r==1,"4-Adjacent","8-Adjacent")))
            )+
  geom_text(aes(label=3*(y+1)+x+1))+
  labs(fill="Pixel Type")+
  theme_bw(20)

cent.pix<-data.frame(x=c(-1, 1,1,-1,-1)/2,
                     y=c(-1,-1,1, 1,-1)/2)
ggplot(morph.2d,aes(x=x,y=y))+
  geom_tile(color="black",size=1,alpha=0.7,
            aes(fill=ifelse(r==0,"Pixel of Interest",
                                   ifelse(r==1,"4-Adjacent","8-Adjacent")))
            )+
  geom_path(data=cent.pix,color="red",size=3)+
  geom_point(data=cent.pix,color="yellow",size=5)+
  geom_text(aes(label=3*(y+1)+x+1))+
  labs(fill="Pixel Type")+
  theme_bw(20)

morph.3d<-expand.grid(x=c(-1,0,1),y=c(-1,0,1),z=c(-1,0,1))
morph.3d$r<-with(morph.3d,sqrt(x^2+y^2+z^2))

cent.pix.3d<-data.frame(x=c(-1, 1,1,-1,-1,-1, 1,1,-1,-1)/2,
                     y=c(-1,-1,1, 1,-1,-1,-1,1, 1,-1)/2,
                     z=c(-1,-1,-1,-1,-1,1, 1,1, 1, 1))
face.pix.3d<-data.frame(x=c(-1, 1, 1,-1,-1)/2,
                        y=c(-1,-1, 1, 1,-1)/2,
                        z=c( 0, 0, 0, 0, 0))
pixel.label<-function(r) ifelse(r==0,"Pixel of Interest",
                                   ifelse(r==1," 6-Adjacent",
                                          ifelse(r==sqrt(2),"18-Adjacent","26-Adjacent")))
ggplot(morph.3d,aes(x=x,y=y))+
  geom_tile(color="black",size=1,alpha=0.7,
            aes(fill=pixel.label(r))
            )+
  geom_path(data=cent.pix.3d,color="yellow",size=3)+
  geom_point(data=cent.pix.3d,color="purple",size=5)+
  geom_path(data=face.pix.3d,color="red",size=3)+
  geom_point(data=face.pix.3d,color="yellow",size=5)+
  geom_text(aes(label=9*(z+1)+3*(y+1)+x+1))+
  facet_grid(~z)+
  labs(fill="Pixel Type")+
  theme_bw(20)

cross.im$thresh<-cross.im$col>1.75

cross.mat<-df.to.im(cross.im,"thresh",inv=T)
cross.df<-im.to.df(cross.mat)
ggplot(cross.df,aes(x=x,y=y))+
  geom_tile(data=subset(cross.df,val==0),aes(fill="original"),color="black",alpha=0.7)+
  labs(fill="Type")+
  theme_bw(20)

cross.dil.mat<-imgStdBinaryDilation(cross.mat)
cross.dil.df<-im.to.df(cross.dil.mat)

ggplot(cross.df,aes(x=x,y=y))+
  geom_tile(data=subset(cross.df,val==0),aes(fill="original"),color="black",alpha=0.7)+
  geom_tile(data=subset(cross.dil.df,val==0),aes(fill="post dilation"),color="black",alpha=0.3)+
  labs(fill="Type")+
  theme_bw(20)

cross.im$thresh<-cross.im$col>0.4
cross.mat<-df.to.im(cross.im,"thresh",inv=T)
cross.df<-im.to.df(cross.mat)
ggplot(cross.df,aes(x=x,y=y))+
  geom_tile(data=subset(cross.df,val==0),aes(fill="original"),color="black",alpha=0.7)+
  labs(fill="Type")+
  theme_bw(20)

cross.ero.mat<-imgStdBinaryErosion(cross.mat)
cross.ero.df<-im.to.df(cross.ero.mat)
ggplot(cross.df,aes(x=x,y=y))+
  geom_tile(data=subset(cross.df,val==0),aes(fill="original"),color="black",alpha=0.7)+
  geom_tile(data=subset(cross.ero.df,val==0),aes(fill="post erosion"),color="black",alpha=0.3)+
  labs(fill="Type")+
  theme_bw(20)

cross.im$thresh<-cross.im$col>0.4
cross.mat<-df.to.im(cross.im,"thresh",inv=T)
cross.df<-im.to.df(cross.mat)
ggplot(cross.df,aes(x=x,y=y))+
  geom_tile(data=subset(cross.df,val==0),aes(fill="original"),color="black",alpha=0.7)+
  labs(fill="Type")+
  theme_bw(20)

cross.ero.mat<-imgStdBinaryOpening(cross.mat)
cross.ero.df<-im.to.df(cross.ero.mat)
ggplot(cross.df,aes(x=x,y=y))+
  geom_tile(data=subset(cross.df,val==0),aes(fill="original"),color="black",alpha=0.7)+
  geom_tile(data=subset(cross.ero.df,val==0),aes(fill="post opening"),color="black",alpha=0.3)+
  labs(fill="Type")+
  theme_bw(20)

cross.im$thresh<-cross.im$col>1.75

cross.mat<-df.to.im(cross.im,"thresh",inv=T)
cross.df<-im.to.df(cross.mat)
ggplot(cross.df,aes(x=x,y=y))+
  geom_tile(data=subset(cross.df,val==0),aes(fill="original"),color="black",alpha=0.7)+
  labs(fill="Type")+
  theme_bw(20)

cross.dil.mat<-imgStdBinaryClosing(cross.mat)
cross.dil.df<-im.to.df(cross.dil.mat)

ggplot(cross.df,aes(x=x,y=y))+
  geom_tile(data=subset(cross.df,val==0),aes(fill="original"),color="black",alpha=0.7)+
  geom_tile(data=subset(cross.dil.df,val==0),aes(fill="post closing"),color="black",alpha=0.3)+
  labs(fill="Type")+
  theme_bw(20)

cortbone.im<-imagedata(t(png::readPNG("ext-figures/bone.png")),"grey")
cortbone.df<-im.to.df(cortbone.im)
#cortbone.df$val<-255*(cortbone.df$val>0.5)
cortbone.close.im<-imgStdBinaryClosing(cortbone.im)
ggplot(subset(cortbone.df,val<1),aes(x=x,y=518-y))+
  geom_raster(data=subset(im.to.df(cortbone.close.im),val<1),
              aes(fill="closing 3x3"),alpha=0.8)+
  geom_raster(aes(fill="before closing"),alpha=0.6)+
  labs(fill="Image",y="y",x="x",title="3x3 Closing")+
  coord_equal()+
  theme_bw(20)

cortbone.close.im<-imgStdBinaryClosing(cortbone.im,dim=7)
ggplot(subset(cortbone.df,val<1),aes(x=x,y=600-y))+
  geom_raster(data=subset(im.to.df(cortbone.close.im),val<1),
              aes(fill="closing 7x7"),alpha=0.8)+
  geom_raster(aes(fill="before closing"),alpha=0.6)+
  labs(fill="Image",y="y",x="x",title="7x7 Closing")+
  coord_equal()+
  theme_bw(20)

full.kernel<-makeBrush(45,"box")
cortbone.close.im<-closing(cortbone.im<1,full.kernel)
ggplot(subset(cortbone.df,val<1),aes(x=x,y=600-y))+
  geom_raster(data=subset(im.to.df(cortbone.close.im),val>0),
              aes(fill="closing 45x45"),alpha=0.8)+
  geom_raster(aes(fill="before closing"),alpha=0.6)+
  labs(fill="Image",y="x",x="y",title="45x45 Closing")+
  coord_equal()+
  theme_bw(20)

get.vf<-function(img) (100*sum(img>0)/
               (nrow(cortbone.im)*ncol(cortbone.im)))
vf.of.closing<-function(dim.size) {
  data.frame(dim.size=dim.size,
             vf.full=get.vf(closing(cortbone.im<1,makeBrush(dim.size,"box" )))
             )
  }
pt.list<-c(seq(1,10),seq(11,40,2))
vf.df<-ldply(pt.list,vf.of.closing,.parallel=T)
ggplot(vf.df,aes(x=dim.size))+
  geom_line(aes(y=vf.full))+
  labs(fill="Image",y="Bone Volume Frac(%)",x="Closing Size",title="Closing Size vs Volume Fraction")+
  theme_bw(15)

disc.kernel<-makeBrush(45,"disc")
cortbone.close.im<-closing(cortbone.im<1,disc.kernel)
ggplot(subset(cortbone.df,val<1),aes(x=x,y=600-y))+
  geom_raster(data=subset(im.to.df(cortbone.close.im),val>0),
              aes(fill="closing 45x45"),alpha=0.8)+
  geom_raster(aes(fill="before closing"),alpha=0.6)+
  labs(fill="Image",y="x",x="y",title="45x45 Closing")+
  coord_equal()+
  theme_bw(20)

get.vf<-function(img) (100*sum(img>0)/
               (nrow(cortbone.im)*ncol(cortbone.im)))
vf.of.closing<-function(dim.size) {
  data.frame(dim.size=dim.size,
             vf.full=get.vf(closing(cortbone.im<1,makeBrush(dim.size,"box" ))),
             vf.disc=get.vf(closing(cortbone.im<1,makeBrush(dim.size,"disc"))),
             vf.diam=get.vf(closing(cortbone.im<1,makeBrush(dim.size,"diamond"))),
             vf.vline=get.vf(closing(cortbone.im<1,makeBrush(dim.size,"line",angle=0))),
             vf.hline=get.vf(closing(cortbone.im<1,makeBrush(dim.size,"line",angle=90)))
             )
  }
pt.list<-c(seq(1,9,2),seq(11,60,4))
vf.df<-ldply(pt.list,vf.of.closing,.parallel=T)
ggplot(vf.df,aes(x=dim.size))+
  geom_line(aes(y=vf.full,color="Full"))+
  geom_line(aes(y=vf.disc,color="Circular"))+
  geom_line(aes(y=vf.diam,color="Diamond-Shape"))+
  geom_line(aes(y=vf.vline,color="Vertical Line"))+
  geom_line(aes(y=vf.hline,color="Horizontal Line"))+
  labs(color="Kernel Used",y="Bone Volume Fraction (%)",x="Closing Size",title="Closing Size vs Volume Fraction")+
  theme_bw(20)

make.circle<-function(nx,r,nscale=1) {
  ny<-nx
  grad.im<-expand.grid(x=c(-nx:nx)/nscale,y=c(-ny:ny)/nscale)
  grad.im$val<-with(grad.im,(x^2+y^2)

calc.eros.cnt<-function(nx,r,nscale=1,im.fun=make.circle) {
  in.cir<-im.fun(nx,r,nscale=nscale)
  cir.im<-df.to.im(in.cir)
  data.frame(nx=nx,r=r,nscale=nscale,
             dvolume=sum(dilate(cir.im>0,makeBrush(3,"diamond"))>0),
             evolume=sum(erode(cir.im>0,makeBrush(3,"diamond"))>0),
             volume=sum(cir.im>0)
             )
}
max.rad<-20
r.list<-seq(0.2,max.rad,length.out=20)
r.table<-ldply(r.list,function(r) calc.eros.cnt(max.rad+1,r,1))
r.table$model.perimeter<-with(r.table,2*pi*r)
r.table$estimated.perimeter<-with(r.table,
                                  ((dvolume-volume)+(volume-evolume))/2)
ggplot(r.table,aes(x=r))+
  geom_line(aes(y=estimated.perimeter,color="Morphological Estimation"))+
  geom_line(aes(y=model.perimeter,color="Model Circumference"))+
  labs(x="Radius",y="Perimeter",color="Source")+
  theme_bw(20)